{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fully Bayesian GPs - Sampling Hyperparamters with NUTS\n",
    "\n",
    "In this notebook, we'll demonstrate how to integrate GPyTorch and NUTS to sample GP hyperparameters and perform GP inference in a fully Bayesian way.\n",
    "\n",
    "The high level overview of sampling in GPyTorch is as follows:\n",
    "\n",
    "1. Define your model as normal, extending ExactGP and defining a forward method.\n",
    "2. For each parameter your model defines, you'll need to register a GPyTorch prior with that parameter, or some function of the parameter. If you use something other than a default closure (e.g., by specifying a parameter or transformed parameter name), you'll need to also specify a setting_closure: see the docs for `gpytorch.Module.register_prior`.\n",
    "3. Define a pyro model that has a sample site for each GP parameter, and then computes a loss. For your convenience, we define a `pyro_sample_from_prior` method on `gpytorch.Module` that does the former operation. For the latter operation, just call `mll.pyro_factor(output, y)` instead of `mll(output, y)` to get your loss.\n",
    "4. Run NUTS (or HMC etc) on the pyro model you just defined to generate samples. Note this can take quite a while or no time at all depending on the priors you've defined.\n",
    "5. Load the samples in to the model, converting the model from a simple GP to a batch GP (see our example notebook on simple batch GPs), where each GP in the batch corresponds to a different hyperparameter sample.\n",
    "6. Pass test data through the batch GP to get predictions for each hyperparameter sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "import pyro\n",
    "from pyro.infer.mcmc import NUTS, MCMC\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Training data is 11 points in [0,1] inclusive regularly spaced\n",
    "train_x = torch.linspace(0, 1, 6)\n",
    "# True function is sin(2*pi*x) with Gaussian noise\n",
    "train_y = torch.sin(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * 0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will use the simplest form of GP model, exact inference\n",
    "class ExactGPModel(gpytorch.models.ExactGP):\n",
    "    def __init__(self, train_x, train_y, likelihood):\n",
    "        super(ExactGPModel, self).__init__(train_x, train_y, likelihood)\n",
    "        self.mean_module = gpytorch.means.ConstantMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.PeriodicKernel())\n",
    "    \n",
    "    def forward(self, x):\n",
    "        mean_x = self.mean_module(x)\n",
    "        covar_x = self.covar_module(x)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running Sampling\n",
    "\n",
    "The next cell is the first piece of code that differs substantially from other work flows. In it, we create the model and likelihood as normal, and then register priors to each of the parameters of the model. Note that we directly can register priors to transformed parameters (e.g., \"lengthscale\") rather than raw ones (e.g., \"raw_lengthscale\"). This is useful, **however** you'll need to specify a prior whose support is fully contained in the domain of the parameter. For example, a lengthscale prior must have support only over the positive reals or a subset thereof."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Sample: 100%|█████████████████████████████████████████| 300/300 [00:17, 17.30it/s, step size=6.50e-01, acc. prob=0.895]\n"
     ]
    }
   ],
   "source": [
    "# this is for running the notebook in our testing framework\n",
    "import os\n",
    "smoke_test = ('CI' in os.environ)\n",
    "num_samples = 2 if smoke_test else 100\n",
    "warmup_steps = 2 if smoke_test else 200\n",
    "\n",
    "\n",
    "from gpytorch.priors import LogNormalPrior, NormalPrior, UniformPrior\n",
    "# Use a positive constraint instead of usual GreaterThan(1e-4) so that LogNormal has support over full range.\n",
    "likelihood = gpytorch.likelihoods.GaussianLikelihood(noise_constraint=gpytorch.constraints.Positive())\n",
    "model = ExactGPModel(train_x, train_y, likelihood)\n",
    "\n",
    "model.mean_module.register_prior(\"mean_prior\", UniformPrior(-1, 1), \"constant\")\n",
    "model.covar_module.base_kernel.register_prior(\"lengthscale_prior\", UniformPrior(0.01, 0.5), \"lengthscale\")\n",
    "model.covar_module.base_kernel.register_prior(\"period_length_prior\", UniformPrior(0.05, 2.5), \"period_length\")\n",
    "model.covar_module.register_prior(\"outputscale_prior\", UniformPrior(1, 2), \"outputscale\")\n",
    "likelihood.register_prior(\"noise_prior\", UniformPrior(0.05, 0.3), \"noise\")\n",
    "\n",
    "mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)\n",
    "\n",
    "def pyro_model(x, y):\n",
    "    model.pyro_sample_from_prior()\n",
    "    output = model(x)\n",
    "    loss = mll.pyro_factor(output, y)\n",
    "    return y\n",
    "\n",
    "nuts_kernel = NUTS(pyro_model, adapt_step_size=True)\n",
    "mcmc_run = MCMC(nuts_kernel, num_samples=num_samples, warmup_steps=warmup_steps, disable_progbar=smoke_test)\n",
    "mcmc_run.run(train_x, train_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading Samples\n",
    "\n",
    "In the next cell, we load the samples generated by NUTS in to the model. This converts `model` from a single GP to a batch of `num_samples` GPs, in this case 100."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.pyro_load_from_samples(mcmc_run.get_samples())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.eval()\n",
    "test_x = torch.linspace(0, 1, 101).unsqueeze(-1)\n",
    "test_y = torch.sin(test_x * (2 * math.pi))\n",
    "expanded_test_x = test_x.unsqueeze(0).repeat(num_samples, 1, 1)\n",
    "output = model(expanded_test_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot Mean Functions\n",
    "\n",
    "In the next cell, we plot the first 25 mean functions on the samep lot. This particular example has a fairly large amount of data for only 1 dimension, so the hyperparameter posterior is quite tight and there is relatively little variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQMAAADGCAYAAADWg+V4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO19d1iVR9r+PXQQUASRInawUhTs3dg1tmg0vZgYTdndZJPdlP3tZ0yy2SRudrMb465fipqiiSV2TayJjYiKiooQwYYCggpIb/P742Z8z4EDHOSIuN/c18XFKe+Zd96Z57mfMk1IKaGhoaFhd6croKGh0TigyUBDQwOAJgMNDY0KaDLQ0NAAoMlAQ0OjApoMNDQ0ANiADIQQLkKIg0KIY0KIk0KIN21RMQ0NjYaFqO88AyGEANBESpkrhHAEsBfAb6WU0baooIaGRsPAob4FSLJJbsVbx4o/PZNJQ+Mug01yBkIIeyHEUQBXAGyTUv5ii3I1NDQaDvX2DABASlkGIEII0QzA90KI7lLKE6bXCCFmA5gNAE2aNIns3LmzLW6toaFRBxw+fDhTStnC0nf1zhlUKVCI/wGQJ6VcUN01UVFR8tChQza9r4aGRu0QQhyWUkZZ+s4WowktKjwCCCFcAYwAcLq+5WpoaDQsbBEm+ANYKoSwB8nlOynlRhuUq6Gh0YCwxWjCcQA9bFAXDQ2NOwibJBA1/ntQUlKClJQUFBYW3umqaNQDLi4uaNWqFRwdHa3+jSYDDTOkpKTAw8MDbdu2BeeTadxtkFLi6tWrSElJQbt27az+nV6boGGGwsJCeHt7ayK4iyGEgLe3d529O00GGlWgieDux630oSYDjUaHlJQUTJo0CcHBwejQoQN++9vfori4GACwZMkSPP/883e4hlXh7u5u8XN7e3tERESgW7duCA8Px4cffojy8vIayzp37hy++eab21HNGqHJQKPeSE1NxZAhQ5CWllbvsqSUmDp1KiZPnoxff/0ViYmJyM3NxRtvvGGDmlpGaWnpbSvb1dUVR48excmTJ7Ft2zZs3rwZb75Z88LeO0UGkFI2+F9kZKTUaJw4depUnX8zd+5caWdnJ+fOnVvv+2/fvl0OGjTI7LPs7GzZvHlzmZeXJ7/44gs5ceJEOXr0aBkSEiLnzZsnpZQyNzdXjhs3ToaFhclu3brJFStWSCmlPHTokBw8eLDs2bOnHDVqlLx8+bKUUsohQ4bI1157TQ4ePFjOmzdPtmnTRpaVlUkppczLy5OtWrWSxcXF8syZM3L06NGyZ8+ecuDAgTI+Pl5KKWVycrLs27evjIqKkn/6059kkyZNLD5P5c+TkpJk8+bNZXl5uTx79qwcOHCg7NGjh+zRo4fct2+flFLKPn36SE9PTxkeHi4//PDDaq+rDZb6EsAhWY1eajLQMENdyMDFxUWCK1TN/lxcXG75/h999JH83e9+V+XziIgIeezYMfnFF19IPz8/mZmZKfPz82W3bt1kTEyMXLVqlXzqqaduXp+VlSWLi4tlv3795JUrV6SUUq5YsUI+8cQTUkqSgSl5TZw4Ue7cufPmdbNmzZJSSjl8+HCZmJgopZQyOjpaDhs2TEop5b333iuXLl0qpZTy448/tpoMpJSyWbNmMi0tTebl5cmCggIppZSJiYlS6cWuXbvk+PHjb15f3XW1oa5koMMEjVtGcnIyHnzwQbi5uQEA3Nzc8NBDD+Hs2bO3XKaU0mLyy/TzkSNHwtvbG66urpg6dSr27t2L0NBQbN++HX/84x+xZ88eNG3aFAkJCThx4gRGjhyJiIgIvP3220hJSblZ5owZM8xef/vttwCAFStWYMaMGcjNzcX+/fsxffp0RERE4JlnnkFqaioAYN++fXjggQcAAI888kidnxHgnI6nn34aoaGhmD59Ok6dOmXxemuvqy/0PAONW4a/vz88PT1RWFgIFxcXFBYWwtPTE35+frdcZrdu3bB69Wqzz3JycnDx4kV06NABhw8frkIWQgiEhITg8OHD2Lx5M1577TWMGjUKU6ZMQbdu3XDgwAGL92rSpMnN1xMnTsRrr72Ga9eu4fDhwxg+fDjy8vLQrFkzHD161OLvbyVjn5ycDHt7e/j6+uLNN99Ey5YtcezYMZSXl8PFxcXib/7+979bdV19oT0DjXohPT0dc+bMQXR0NObMmVPvJOI999yD/Px8LFu2DABQVlaG3//+93j88cdveiDbtm3DtWvXUFBQgLVr12LAgAG4fPky3Nzc8PDDD+Pll1/GkSNH0KlTJ2RkZNwkg5KSEpw8edLifd3d3dG7d2/89re/xYQJE2Bvbw9PT0+0a9cOK1euBECLfuzYMQDAgAEDsGLFCgDA119/bdWzZWRkYM6cOXj++echhEB2djb8/f1hZ2eHL7/8EmVlZQAADw8P3Lhx4+bvqrvO5qgufridfzpn0HhxKwlEW+PChQtywoQJsmPHjrJ9+/by+eefl4WFhVJKKb/44gs5ffp0OW7cOLME4tatW2VoaKgMDw+XUVFRMiYmRkopZWxsrBw0aJAMCwuTXbt2lYsXL5ZSMmegrlFYuXKlBCB3795987Pk5GQ5evRoGRYWJrt06SLffPPNm5+rBOK7775bbc7Azs5OhoeHy65du8qwsDD5wQcf3ExUJiYmytDQUNmnTx/56quv3iyjuLhYDh8+XIaFhckPP/yw2utqQ11zBjbfz8Aa6P0MGi/i4+PRpUuXO10NDRvAUl/e1v0MNDQ0/jugyUBDQwOAJgMNDY0KaDLQ0NAAoMlAQ0OjArbYEDVICLFLCBFfcbzab21RMQ0NjYaFLTyDUgC/l1J2AdAXwHNCiK42KFfj/yjeeecddOvWDWFhYYiIiMAvv9zeM3mGDh2Kugx17969GxMmTLD4uRACn3322c3PYmNjIYTAggXVnhzQaGCLDVFTAaRWvL4hhIgHEAjg9kyg1vivxoEDB7Bx40YcOXIEzs7OyMzMvLmXwd2A0NBQfPvtt5g1axYArnMIDw+/w7WyDjbNGQgh2oI7Jevj1TRuCampqfDx8YGzszMAwMfHBwEBAQCA+fPno1evXujevTtmz559c8HP0KFD8eKLL2Lw4MHo0qULYmJiMHXqVAQHB+NPf/oTAO4R0LlzZzz22GMICwvDtGnTkJ+fX+X+P/74I/r164eePXti+vTpyM3lMaJbt25F586dMXDgQKxZs6ba+rdu3RqFhYVIT0+HlBJbt27F2LFjb36flJSEMWPGIDIyEoMGDcLp0zxiZMOGDejTpw969OiBESNGID09HQAwb948PPnkkxg6dCjat2+Pf/7znwCAvLw8jB8/HuHh4ejevfvNRVb1gc3IQAjhDmA1gN9JKXMsfD9bCHFICHEoIyPDVrfVuI0oLgYSE237V5uRHzVqFC5evIiQkBA8++yz+Omnn25+9/zzzyMmJgYnTpxAQUEBNm40judwcnLCzz//jDlz5mDSpElYuHAhTpw4gSVLluDq1asAgISEBMyePRvHjx+Hp6cnPvnkE7N7Z2Zm4u2338b27dtx5MgRREVF4cMPP0RhYSGefvppbNiwAXv27Kl1/cW0adOwcuVK7N+/Hz179rxJbAAwe/Zs/Otf/8Lhw4exYMECPPvsswCAgQMHIjo6GrGxsZg5cybef//9m785ffo0fvjhBxw8eBBvvvkmSkpKsHXrVgQEBODYsWM4ceIExowZU3PDWgFbHbzqCBLB11JKi7QppVwspYySUka1aGHxqDcNDbi7u+Pw4cNYvHgxWrRogRkzZmDJkiUAgF27dqFPnz4IDQ3Fzp07zRYdTZw4EQDd9G7dusHf3x/Ozs5o3749Ll68CAAICgrCgAEDAAAPP/ww9u7da3bv6OhonDp1CgMGDEBERASWLl2K8+fP4/Tp02jXrh2Cg4MhhMDDDz9c4zPcf//9WLlyJZYvX35zmTOAGpdEp6SkYPTo0QgNDcUHH3xg9mzjx4+Hs7MzfHx84Ovri/T0dItLtuuLeucMBNdxfgYgXkr5Yb1rpNFo4OQEhIQ0/H3t7e0xdOhQDB06FKGhoVi6dClmzpyJZ599FocOHUJQUBDmzZtntvuvsr52dnZmltjOzu7mtmaWlj6bQkqJkSNHYvny5WafHz16tE7Llf38/ODo6Iht27bho48+wv79+wEA5eXl1S6JfuGFF/DSSy9h4sSJ2L17N+bNm1fl2VTblJaWWlyy/ec//9nqOlqCLTyDAQAeATBcCHG04m+cDcrV+D+IhIQE/PrrrzffHz16FG3atLmp+D4+PsjNzcWqVavqXPaFCxduLmdevnw5Bg4caPZ93759sW/fPpw5cwYAkJ+fj8TERHTu3Blnz55FUlLSzd/Whvnz5+O9996Dvb39zc9qWhKdnZ2NwMBAAMDSpUtrLd/Sku36whajCXsB6L21NWyC3NxcvPDCC8jKyoKDgwM6duyIxYsXo1mzZjd3+2nbti169epV57K7dOmCpUuX4plnnkFwcDDmzp1r9n2LFi2wZMkSPPDAAygqKgIAvP322wgJCcHixYsxfvx4+Pj4YODAgThx4kSN9+rfv7/Fz7/++mvMnTsXb7/9NkpKSjBz5kyEh4dj3rx5mD59OgIDA9G3b99ad4uKi4vDK6+8Ajs7Ozg6OmLRokV1aAnL0EuYNczw37qE+dy5c5gwYUKtSvzfBL2EWUND45agyUDj/wTatm37f8oruBVoMtDQ0ACgyUDDAu5EHknDtriVPtRkoGEGFxcXXL16VRPCXQxZcSR7XbdU1+cmaJihVatWSElJgZ4yfnfDxcUFrVq1qtNvNBlomMHR0RHt2rW709XQuAPQYYKGhgYATQYaGhoV0GSgoaEBQJOBhoZGBTQZaGhoANBkoKGhUQFNBhoaGgA0GWhoaFRAk4GGhgYATQYaGhoV0GSgoaEBwHZbpX8uhLgihNC7R2ho3KWwlWewBED9T3HQ0NC4Y7AJGUgpfwZwzRZlaWho3Bk02BJmIcRsALMBnkdnK0gJ/PQTkJ1tfBYeDrRta7NboLwcsLOSNjMygOho1ktKQJ29YekMDvW92kdESsDZGRg5EjDZbr9R4MgRICXFqGu7dkBY2J2pS24usHcvj2ozbVfT9jRtd9P2Bdifo0cDrq4NW+/akJwMnDxpyI6XFzBwoPWyV1/YbKv0ikNXN0opu9d2ra22Sr9+HVi7Frj3XsDHh59JCezbBxQUUKluFWVlwIYNfF1aytOF7r3XslIDwMWLFNDWrYE+fQCHW6TZrCzgxx+Bnj2Bjh1vrQxborgYWLUKiIoyP10pPh44cQK4776GE9acHGDTJvb1wIG3rszFxcDWrYCfH9C7t23reCuQEti4kc/Vrx/r98MPwJUrwNGjQP/+wJQpQB03LrKImrZKh5TSJn8A2gI4Yc21kZGRsr7IypJy2TIpy8stf3/hgpTr1t162Z9/LuWNG8Zn165JuWSJ5fv9/LOUu3bd2r2qw9dfS1lWZtsy64rycik/+0zKwkLL32dns00aAidOSLl2rZSlpbYrc8MG9vWdxnffUb6kNOQsP5/vS0qk/PRT27UzgEOyGr28a4cWt24FHnqoeksdFAR4egK1HJhbBVeukKUfewxwdzc+9/ICRo0Cfv7Z+Cw1NRW9ew/BtWtpGDq09rKLi2ndrl8H8vLorlaHceOAzZvrVndbY/t2ekMmR/2ZwdMTGDoUqDhK8Lbhxg3g9Glg0iTbhk/jxgFbttiuvFtBQgLQvj3l6/JlegSPPmp4PQ4OfJ+VBRw+fHvrYpOcgRBiOYChAHyEECkA/kdK+ZktyraEc+eo7LW5p0OGAN99B8yYYV25V64AO3cCDz5YlWSuXAH27GGHpKXRjf/b397CoUN78cMP8zFp0icoL+ex47/+asSnpjGskxNdPTs7oKQEyM+3TAju7kCPHlTCjAzA1odWX70K7NrF+hQUMMfSubP5NQUFjM1ru3ebNkaOpA5nk1oNKYHvvweqO/j42jXmjIRgm7ZvD0RGWle2nR3QrRsQFweEhtquztZCSsrTgw8yH7NzJ9C0KZ8XAJo0Yajr6MiQMT6e9bVFuGAJd+XxasuXAyYnXdeIw4dpwYKDa77u0iXmGqZPryrUu3bRqo8aRcveooUrSkoKq5Th6OiCZcsK4OLCnMONG1SqwkL+qYSXgwMVvlkz1s3e3pw8fH3pOWRnAzExtL79+5t7KreKgweB9HRgwgTjOWNiqFSjRxvXrVpFS+zoWHuZ6enAqVPAsGH1r19l7NkDdOgABARU/e7UKVrWyZONZzlzhrmMyZOtv8fatXW73lbYtYvKXVwMrF/PnMGUKUabX79O72z6dOatvvuOMjOuHsca15QzuOs2RL1woW6JtchIYM2amskgORk4fhy4//6q323fTobOzwcWLKC7FhmZjEuXXkZKylpImQ87Ozd4eU1BQMAC/P3vgJsbSUApv50dCcDenu9LS0kWRUX87+hIUggIAAIDSRze3rRWUvL63btp+ZycSA5NmtS15ahYHh50/U3Rqxdw6BBd8c6dWScHB+uIAABatmTZtwOZmcCgQVU/j4lhe0yZYv55x47sq9hYeleNGVlZVPgffwQ6daKxMYWXFxARQc9nyBB6Ddev37763HVkcPhw3Vnczq764cFDh9jApmWWl/Pz5cvpCXToQNb28GCIEhLij4sXPSFlIYRwQXl5IRwdPREU5Ad3dyqTuzsVukULYNas6t3t3Fzg7Flaubg4lp+fT3Jau5bXZGby/+DBJBoV57Zpwyy/Ne75kSMkkogIy99HRfF5g4OZAxgwoPYyTdGmDUdUgoLq9ruakJ9vuP+mxJSYyHarzhMJC6MB6NCB1x08aD7UWF5Or0d91qYN292Ww9G1ITOToWdJCQmsMhEoBAezXS9fphFYtoxhaHAwvUcHh+pzOnXFXUcGQN1jU2X5TIeRpGSCLjCQcVlxMZUgJYWK6epKy+zrSxcuPp65guvXeW1GRjq8vefgmWdmY+PGxXBwSEWPHkaOoKSEFv70abp1JSW8r3J5XVzYiUKQDHx8KMT33kvX0d0dSE0FvvySVre0FPjb3wB/f1qI7t2pLGfO8L2KLS0hJYXCN3w48OmnwFNPWb5uwgS2SWlp3fMUkZFsJ0UG335rfa6mOuzezTbfs4d1B/gcJ09W9QhMUV5Oq/rmm8C0aeaKD9AiL10KzJzJfggPZ9LYVmSQlsZhy+qwejWN2vjxJKGxY2sub9gwhm3Tp7PcuDiSwbFjNA7VEXxdcVeRQVISE0R1RWAgrYNCfj4bt08flnnuHDvQxYXXDh0KPP888MsvZOTSUrrlPXoAEycC//oXUFS0BunpHF0oLl6Itm3pEdjZUfCaNuWYcUgIiUUlDVevJrFcu8b7DRvGz52dSTTKG2nalB3drh1wzz3Ab35DRfj8cxLb+vW8pmlTktbly1Tg0aONBFNqaipmzJiJ++//Fs8/74ddu2rOxnt4MHQx9aCsTQza2Rl5D4CezdSp1ocallBYSKLMyCAZKAJ/5BHL1+fmMqwrL2doUVDAtqlc/2bNOBK1ciWTd8pztBXefpshpaVEX34+rfu99/LZrMl9CcG+vXLFaOeSEsrQjRu2q/ddRQbHj9cv0SMllff77yksqakUjEuX6BaXlAAffkiBcnentR45kkRw6hSve/11EkaLFuwIV1fgtdc4GaZ3b/6ub18qd2IiLbfpjLIXXqCSA8D588CKFSQkgPF6RATjx19+IUHt38+JJ99/Ty+lY0daCDs7WrPEROY8mjcnUZ44QQ9jzBjgrbfewt69exESMh/AJ7h+nfVW+YrERHoNShGEYBsXFdH9rKxE6jkAtsnQoeaTqzw8SGSennzGS5du3dpmZpLoHByo1ABj6zFjqtaroIDt7+5Oa6sIaORItpElL8LRkUSdlERvrfLMxVvF8eP0vPbsqTrp7bvvaJRKS9lnTz1l/VDpkCE0YP3705js28f6FhfXr76muKvIAKh7Z0lJS3rhAt1GHx/gH/+gsl24QPc2JYVjuampFKjevemGX73KeK1lS1qYrCx6B66uVFAPDyrlli283s2N3kB0NEcCXF1pqVVMd/06hUQJt4cHMGcOCUpKEtSnn/L3zZpRYAMDqZze3iSHn36iV+DsDHTpAowYQaXbv5/us709cPmyK6Q0Rjs++2wRPvtsEeztXfCXvxTg2DG6xt27kzhMPYHSUv7VNNsSoNJv2MDfjh1Lj8bLi883fjzj8AsXbp0MoqNpWbt2JRGnpDDn4etrfl1MDL+bOJHfm0KRgso5VFb2qCiGMx06sC1OnLi1Icblyyk3e/aw/DFj6KWY4vBhjnwsWAA8+yzDx7p4uUJQJtzcWLaUJEpbDgbeNUOLv/5KtzEkhEMyY2pZI5mdzXFbKZkzcHQkCcyYQYsYHg588QXw1Ve0lJ07UyAuXqSXcOUKBc/Xl0LWty+TcFOnUsFfeokuf04OXf24OFrn+fOpnIcOkbWPHWP5c+dWVYwbN3jdjh0kHnd3hgSdOzP2VlYjO5vZcTWxqbycgrtkCdulWTN6C61b00u5dCkVFy++jKystQDy4ejohoCAKRg3bgE6dvTDsWNUgJISClnr1vQYSktJcgMH0vIPHGjUtXISTyE/H/jjH+nxtGgBvPEGwyg1YaryyIW1WLeOz+noyPa8cQN4+WXzazZvZv6lppg5N9fIOXz1FZO5poiLI6F17cp71tXzzMxkm3Xrxr785Rcaj7g4ErXyQJYsAd5/n+09ZQoJs0OHug3Hmk6RN02ITpxofRn/FUOLauz4yBFaHFPs2EElAqg4P/xAF3PiRCpUaSnw0UcU6F9/pbAmJNC69u5N17tnT8ZybduSfV95heUq61lWRmH54x9JDqNH020tLASeeQb4/e8Z2y5ZwiSfck1nzOA1777LzwcMMKyPhweFITaWyvfaa0ZS8IcfeE+AhPbLL+bxoZQkh3vuYey4fDnzBr6+wKFD/ujb1xOxsYUAXFBSUoiyMk9ERPhh9mw+x6RJLKe0lB7F2bN0P5s0YRt8/jmfp2VLejuLFjFvURlxcSS6EyfYlgMHsl3qi6Ii3ByZKSur6k5v3EgSqO1sUXd3tumpU+wHU+zbx/5Ys4bKXJvXeeZM1WHtvXvZdmrE5+pVysK6dezfIUOA555jP6kwzM+P/aSS0ZU9mupgb8+2UOv8du+uOlmsPrhryABgY6akmGdqc3KAbdvI/Fu30ppMm2a4vjk57Bw3N3oUy5ZRKMaOpSKOGEGB/vlnlr1lC3DgAP9mzjTK2bOHIcahQ8DChSSTDh044SYvj0L144+sw3ffmdfbxYUhypYtzBMkJDCuHjaMnkh+vrEC0M2Nry2tCLRkZcvKKIje3iSEuDgK9oUL6QgOnoPi4tm4fn0xsrNT8dVXJJXx443fOzjQSwoP5/vNmw1SaNWKlnflSraxJTK4dIneUteuwNdfM4zatInW8FYWMCUk8BmOHaOXsWMHlc3USzl4kAph7SHDnTtTcUyHPYuL6dkNGGAoWW345z/p5psqb3k5CaF/f5LnP/7BNpw2DXjnHXomo0aRRH76iSQ2ejRlecQIflaXBXX9+jHfcP06jYAtcVesTbh0ibEzwEY0ZfBz56hYX35JyzR+PIUwOxuYN4/TWE3nDAwdSpd2yhS6tcuWsUOHDePowvLlVPrJk82F+do1xnkBAYyH4+MpkC1b0rJPmEDFlrIqGSiMHcu69+nDTl29GvjkEypkaWnNbVB5qe7BgxTM//kf4OOPabVatuTz9+0L3HPPGhQWLsT16+F46KGFCAtbg6IiCu6HH9KDMUV5Oe8xfDjJ6//9PxLAnj2c3+DlxTY2FUAVSinMnMkYPiCAymwJq1bV/JwJCWzf554jMZ49yxEVOzu20aVLzN3UZfl0t27sG1NcvGgMn/bvz5yL6XJnS+jXjwZF4dQpKnliIsnznXfobdrZkTAyMii37u4MVYuL6d2pGaeenpRJha++qj0HEBBAA2TqNdoKd4VnEBNjHheZNsDPP7OBH3mEibktW5isWbeO7y9coJL8+c9k6O3b+dnOnRT8V15hp7z6KjtpzBgKvilycqggiYkcCjKNn8PCqJhpacZw0d/+Rjfc0mSQ8eM5gjBuHBXw/fdJQEePGjMUr1xh3du1M571wAFmvtPT+b5zZ4YggYEGaa1fz+9dXFjn/Hwqz5YtDFFcXalQmZnAX/9KL+F3v+PoRUwME2pHjhiTfbp0Icn97//SGxs/nuTQsSMV7MABfq9gb08r9/nnLMOSoP70E9u4pqnVqalUHoChz4ABDOOOHGHfWrMoTEGtKcnIoAencO4c27e8nKSwdy89jdRUy1OfAdbZNDGoZmw2aUKiUKHMPffQo2ndmn3y6qv8PC2N3uTly+wPgLKp5iXExlIumjev+Zl8fFh/NVnKVutC7grPQEoKvBrHV5/t2kWh9/WlpVu1CnjvPca3Fy6QJCZPprJGRdFaXbrEckaOZKd5eNCylpcDL75oKJspfvqJoUXnzuzwXbvottrbs2NUyPD00+xoO7uaV0vefz8Jo2dPXuvpSWs+bhxJz9+fhDRmDIWlrIwCN2gQcw/vvstVlaaLtZTixcTQKsXGMnT49FO6pU2b0qo7ObHO5eVU5tmz+fzJyYb3FRTE9hsyhIrUujWJ7t//pneTnc14u7y8aiwfEEClCwggMZoSQkoKlTw1tfb+BngfNzeSqrc3+8HFhSRmjVufkECLP2UKvcYVK4yhuLw8hjYq/yQEQ5vTp2suMySE5QL06NavZzmmnp27O2XtwQeNBG9xMcn8hRfY1oWFbD/llaSmUh5qaxuAcpCTw3tmZ1cdubhVNHoySEoyYsPERHaGEMwsnzrFBunQAVi8mFY+Lo6M//nnwFtvUcH++lda8vJyEkf//ozpfvyRrq+LC5W7dWsKbGWUltJy9O5N8lBWz8ODnXr0qGGBk5J4r48/pqBs3MhhsrNn+TtVXuvWfJ6CAlrsFi2YyPvpJwqHnR3zDK1akdDmzKmaADPFuXMsq0kTekXvvkvy6tWLG5D4+7O8oCASVrNmxrDlJ5+QTHfu5H2jomiBmzdnQqysjGTVvTvwzTf0GGivZLsAACAASURBVIqLqxfcxx8nGRUU8PcKysO7fNm6vt++nQri58f2TU8nKUybxratCUlJVPTJk41Vpj160AtT6NzZUP4+fShP1SmWIqjQUPaTCquio417mMLdnR5Ykyb0TrZvp1wKQc/m1Cm2n5pEpMqxhgzs7Nh3KSnsR2vbs9ZybVPM7UF5ORspIoIK8tVXVHalmCUlTG7NmUO2DgpiZnjdOirf1atUBnt7NvL999N9v3CBli0+nh2gsrNubobCKmRlsWMvXaKF6taNyr9yJQU+LY35Bikp6CNG8P3p03w/dizHk3Nzec8NG5jL8PKihXj/fWbiT55kueHhfLYHH6QFUKGGt7e5YlVGfDzbIz6eowFqGa9a+TZmDHMMs2Zx1OLqVV4bGEjySksD/vAHXpuebngcoaFMejo6UpAfe4xJ0gsXWFdLa+zd3emBuLnxGU3RtKl5nGwKKdnePj58nZ1NkvT3Z5/27m2s+2jevHolUHKjEnMnT7LfHnyQ9T1/np87ORmegr9/zd5cYaExo9DenmTj70/C69rVfAbjyZOUgY0bWYd9+/g3ZQqfycuLQ8ErVhhtfO4ciaMmwjdFr17UhdJSjo7ZwjtotGRQUsIklhAcZmvfnqRw8CCTP/PnU7ivXKEbPGECw4OMDLpo5eW0qikpvLaoiGWlpzP2DApimLB1K93h6rBvH+vi5cU4cM0aEkqfPiQSb29jg4+EBBJOYSHrBVBwfH3Z4SNHMv4dMoQZ+N/9joI9ciQTnb6+ZHqFPn1IFNYgM5PC4epqTGZJSKDQ9u/Pek2YQOX08uJz9+9PAiospOV0dORMx7feYkLyyhVaT1NvSQjmTU6cYO4gK4vEa4rsbD6vGsEBjEVMNcW2167xr2dPehFt2pAYmjdnXQYMMMilprZZv948x1RQwLoEB7N9qpviUlPdsrNJZADlat8+1rFVK8OYKJw+TfIJCeFQ4smTJDWVMHR2JvEcP86RKeU5ArUnBIuLSTJq6PK551iuLZa3N1oySEtLxebNQxAWloayMg7ZvP46XfvDh9k5M2YwNOjVi5bMw4PKlJFBCzJ9OhVPTa1dvpy/HTKEHdivHwVKWT1LUHsS5uXxen9/ups5ORTckBAjl3HkCBOKrq4UYEvLerdtM1+hFhpKwVLDW6arKz09zTd6rQmnTjGZ+vHHFP4vvqCyDh5sPNurrzKU8fBgLmLUKOAvf2H7RUez3t27k9hWreIc+0uX6KafOWPcq7ycYdXOnfRAoqPNPaqcHCqyiwv7AmCcX9umI0lJJCpXV97X0ZFksGcP6xUQYIxmuLhYtqJqVMXDw/jMdK+Itm1JnKbJRAUvr+rn+puSQUQEifbCBbaNk5NxD0U8arTgD38wn5q8ezfJKCqK4VuLFuzvjh3pbVZHcOXlHLLdvJlGLCqK+Z5PPzV/xvqg0ZLBm2++haNH92LcuPlYuJCKXFTEBp4xg8NY3buzgzMyqKQJCRTs5s2NmWY3bhhLge3taemGDmXyLCvLmEOv8hJKccrLaSXj4ylcppuBCMFyExLYKQAJQsWfI0cyPvzoI/NnSk6mMJoST1ERvYT169mh589TkUzrUhvUXHcvL1qOc+dombp3N8pRMW7btkYStGlTEuaXX7I9i4vZFqGhbMPVq5l7sLc3t6anTrHsmTNpoSZONB8yzM5muOTsTIW5eJEkXdvzpKXxOkWuap1DRgbrevly1Ux7aSkJbvVq9vOiRfRkduywvPioVSuWbZoEVCtK+/Th81vaM0DVBWB7JCfTA6tskdXeA2r3JB8fejoJCQwL3NzoWSliGT6cZNG3Lz0XPz/mZUyHQmNj2b6DBtEbdXUlKQ4eTI+ibVvrjUZNaHRk4OrqCiEEPvtsEYByXLy4CNu2CaSkuN6M+7p3ZydmZtKazJxJlo2KMjaAOHaMArJ/P4fEpk7lb/LyqKhXrtCiqeyw6c6/Z86w44Sgkjs4MEnp5WUIdH6+ufUpLzcsVf/+FGw3N2aTASNJ1KdP1Wf28aFgJSQYSVJLcHe3bLk2bOB9pk+n0vbpQxe2Xz9+37694Ya+/jrw2Wesb2QkvRkhOF8hOJhtlZFB4e/alYq2di3bOTaWZajVoypkWLWK1krNPCwupsA6ODA0WLjQmCFaG4Qw9lNQ24Ipwo2PJ1GqzH1ODvNIgwax/9VzTJ/OfMaSJbTUledodOnCZ5SSba36yMWFxB8fX7Vepp6BesZr14wRGIX8fCNX0qEDSaF1a+a5RoygzLi5URb272dye+JEelnffUfZc3Xl9evX0wg6OzPMU2QEsA0cHEgiSUnm4eWtwlbHq40RQiQIIc4IIV6tT1nJycl48MEH4erqVlG2G1q1egiPPXYWAwaws5o1I9va23PRR2AgBeOzzygcZ87QA7jvPuYT1JwAFf8KYUxT7tqVCtGkiTHt98oVko4a0lQ70tjbG/Pgjx83BDw9nULv40NlFoIdGBpKRVI7FYWFVe/OhYZScM+fr36Puy5dqgpqWRmFpqyM8wV8fGgprlwx1kIEB7NeAC2P2nYtN9fcAk2dys9Hj6ZCXbrEtnZwIDn++KOh8ErB1Bbyv/zCdjYtr00blnflSu1ewdWrVDgfH75WW9+fP0/yBqhogwZRkc6coYeiNq49cYI5F6UwPj5cfPbFF1XXVOTmkuR37GDZpknOpk2NqcWmUKM+AMOxTp2MVa8AZSMtzai3wpkzlJ9hw/h9djZl1N+fRDR9Ogl48mTKoMpFFRczLEhI4H2rG0odPZrX2GJEod5kIISwB7AQwFgAXQE8IIToeqvl+fv7w9PTE4WFnFcPFCI01BP33++H++9njFVQQCY9dIgd+csv9BYiItiZr7xieXOOrl0NYfX3J4F06sROSk1lJ/XtSxe7rMyYMabumZFhuPBqDwKArl7//gw/tm6lwvfta7i6avZjdHRVF9RUSSIiKPwq2135ez8/1rWszHALly7lM4WGsg65uRQe0zkATk6GKwzQ6h49SmUSwlzQXn2VFs/ZmZONsrLo8URGUuBee41KZBpzN2/ONnN0pJus3PMePdhegYG1Z7s3bqRCd+liLIu+fr3qKsWWLVnvU6fYxkIwPEtLoydUXm4os4MDrXFsrDkJFxUBTz5Jj8re3jyc8POzPM3XdGLPjh3sK5WUBihPmzaxbXNy6AFIyZBt1Cheu20bSfmRR9gu06YZnqJaKLdtGz2EyEjgP/9hYvk3v+HsRkthz+nT7O/qJkrVBbbwDHoDOCOlTJZSFgNYAWBSfQpMT0/HrFlzMHBgNPr2nYOkpLSbGdcOHYzlu5Mm0cW65x4KjdpfsKYNNdSa/JQUumtC0NrExbGTFGHExLAsF5eqqw2PHjWsFWC4zFFRFILt22mN09KoNL/+SlcxP792Bu/Th8Kr9jgwFWIleAkJJK6vvuLoSnY2E6wK7u41D1HNmkViKivjs+/ebf79vHmM9X/9lUOfL71E4u3Vy9izb8ECI2wASKolJSQllV9Qm7xMmsS6Kjg7V61ffj49giNHjD0P4+PN1yQAfNZr18w3El27lqT13nskq3ffNRTn2jW62OvW8X3z5rzWzo6hxNatVZc1Vx4OrVxPKamATzxh5DCCgmgs7O3p/qvJRE5OVOyMDNYlMtK4n6nbD5AMTp6kh6DyPR07UsbbtmU/LFhg1GHjRmPVaV2PBLAEW0xHDgRgOvM7BUCVyLgux6utWbMGN26o9foL0b49BaNnTyOGbNeu6sSTgoLap3K2b28IcVgYhy9nzTIUPjCQybXCQoMsNm2iK6c6MSnJMhMLwbyCvT1zFhcv8rOXXqJieXuTDFRHA1WVXUq6iT//TII6eNBYtScEvaCYGFrafv2YEDt/3tw9bdPGyOJbgqMjf1daSmLJyjL/3s6OWfBFi/i+XTsK9uTJVLJmzZhcvXyZBPCf/7DMYcPoBTk5GWsKRo+mtVPkBhgejinJNm1K66l2eSopMXaNSkriiIsQ9AiHDeMks0uXOKls5kx6JLNmkah++YWzJcPCSB5RUSTl2FgqjhruHD+eeSXTWZTVrRNRfb9zp7Hq1M7OMAqmnoMKKd56iwp7+jSV2tHRCNdMISVzM56eDE/z8rjPwsSJxoQ0OzsavbQ0hhaJiTRcTZoAV6+mYvLkmVi79lv41bTfWi2whWdgKRqsEhlLKRdLKaOklFEtrNhgz8ODyZX8fApGRASt//jxdBVbt64ah546ZXknXVP4+TE5U15OgomIMFckZ2danKIiComjI129zEw2vunkE0sYNIgd98ADFNSgIArq1Km0ekVF1f/Wy8uYs969O62zoyOFq0ULkonyfkJC6JKePWuMqpSUsM6urlWPHqvcVoMGsb3S0kiQSUnm1yiLDtBSHjxIpW7TxsilxMby2f78Z0PIp09nW6k2dndnmzs6Ghn8gICqM+0yM1lntevzjh203Js2GRvWenpyFuXw4Sz73nt5TUwM8wM+PrzfiBH0zG7cYKbe3p7XJSaahx0tW7IupudcALzeNKwCjO8VYQUGmu91GB1N70iN2sTFUXZTU9nOI0ey/yrnIzIyOJoTFUUvaPJk9sXUqfRa1q0juS1fTgPg5UUDtmABPaL//V/A3v4tHDy4F/Pnz68qVHWALcggBYDpnritANQ7nVFWRgvTowcF98ABKodyL11dq+4fV1RU+06xJSWGGz9tGhu/soehVsY5OdF6tWjBTu3YkValpg0pvL3pDrZsacyku3qVAtKxY9UJOqbo1YvPm5nJiVYvvMBn7NGDv+3UiULTuzct+5IlVDLlSp85Y2zl9dxz5mVXTlxOmsS65OXR5T5ypOo1XbrwPmqtwtatvNdbb5G0Cgt5/127gL//nQKrpjPb29MqAvRgOnUyxsRVG5ni4kX23+DBrMeNG1Tejh1ZXlkZDUFmJq3oX/7CMtautXzWxciRxojI11+zzHvv5W8r76HYvTstvkJEhOU5ImfP0ghduMBrFCkCxnBobCy/W7mSCexFi4y9LYYONUJHwMh7PfKIQVLe3nx2R0f+btQoyuLYsRxSnzaNuRIpARcXV7i7C2zZsghSlmPRokUQQsD1Fg+htAUZxAAIFkK0E0I4AZgJoJaZ47VDCLraYWGMN3182PH/+pchtKbDOgUF9CbUlmIKhw7RegG0INHRtCyqE+3t2bimx6ZdvMjOKSwkMQQFGUNUalJJTWjVylgQFRrKe23aRAK6dMm4rrLyOTuzPps30/UVwliyqjyDoCAqQU4Ohcbb23iW5GS6rdeuVV3i27Sp+Vi0oyNd4sOH+exqzUVljBjB+/ftyxGU1FQq2rvvGuFObi6Tf7t2cWivsJDfpaYa5HfffSQr1Y6mz15czHv4+NDr2r+f9VHTq7/8kuR4+DAJ4v77WX81alQdJk7k78eOZTjg5sYQx3Q/AkU0586RjEpLSYJff02FVhCChOnpyb4w3fsxN9fwaC5cIBm3bctnMt2Q1cuLcrV6NWW6UydjWbvpfTp0YBtIyVGcsWNZbzs7Pv/q1fTU9u1LRufOD0IICqSbmxseeughnK3J4tSAepOBlLIUwPMAfgAQD+A7KeXJ+pZbXs5Yq1MnWru+fdnxbm6czFN5ddnevZx/b7oL0rJlxhqAb75hRzo4sNO/+sqYl96+PT9TilpcbIQCZ89SIaXk99ZkbXv3pmsHUHDT0/m+bVv+HT3K7/LzzQ9DSUigC+jgQCt84waVIi6ORLh+PRW6uJixsZub+bRbSzsCKbRpYz5KsWIF2/TiRWOS0PHjln87bRoVISCAirR0KS3uU08x4ennx7oeP848wpo19ApSU+ndlJTQfZ87lyM9lXHhAvtbhWsHDvD5mzfnZC+17qN5c5LJ7t1UiuHDjTkCANvZNAxwdTWGKnv1YpsNHWoYB8CYbNa5M/MMq1bRg3F0NA8HleVXaxJMsWePEZ6WlrIeV68ay5sViovZPnFx3J258hwFhZYtmeBev55eRkCAMfGooIDvt2wBJk70x4ULPL/DyckFhYWF8PT0vOW8gU3mGUgpN0spQ6SUHaSU79iiTAcHMnS7dnSjdu/mhJn4eHZGbq75oZm5uWx8RQaFheyII0colA88QFKJi2NyrEcPc+YfOZKW7dw5Wl1vb7rPCQlsfF9fsrGlSUOW0KoVXVoXFyY+T55kfNiqlRFfq/H0+HgOJ50+TVf2wQcptNHRdL0PHODn995LRUtIoGLcuGG5PpbG9AMDDbIrLSXJPPooSTA1lUJWXm45eebszNAoKYkKWFhIhSkoYD0jI9mmJSVs7/R0emRFRfSGFixgOaNGsX8qL1RKSjL67x//YB/88Y8kq++/Z91LSxnaTJ1KV/2jjwwXWsHSmYnBwVSkNm1oDOLjjc1PFZo1o9eYn8+Qw9WV/W/qvaiRjfR09qHpcnqVMMzLY9/a27Nv1DFoZWWU1Q0b+AwdOphv3W8JKnTbsIFlurnR43vqKYYLaqt9B4d0+PvPwcqV0ZgzZw7S6jGs0OhmIJpCZf7btGGDZmQwo5qSwsbv3JnsWVrKDlDHmqWlcYuqggK6pyrBprblUvF75eEttWlmSAhDlOBgCq+zMzuwtNT6KcL9+hneS3AwiUnFoV278rvLl+kO5ufT8hUUGMrt6kqCkpL/t26lYBQUUDHz8qoukFF1szSxydHRSIolJtLjcnMjsRYV0VqGhzNxZwmBgVQae3vjqK8NG9g3QUH8/9prLOPsWRKPgwPbeO1a1h/gNS+8YF72zz+zr3ftYt++9BKThsnJtKbNmpEYFBISaBi+/JL1UOsmKntaCr16Ufn69iUJtGtnTga+vuyPp59mFj8sjIpmmuRUKyiTkkh+CQlsw8JCI0917Rrvr1baqslQX39NebjvPuNsjkOHqp+AlpHBERp3d3pz585x8tQ777CPXV3ZJoMGAVu3rsHFiwsxYUI4Fi5ciDVr1lgu1Ao0ajJo2dJwbSdO5DiyoyMZ84cfqJxNm3KvPB8fWtHNmylUHh6MB0+eNKYmp6cbU1t79mQjm2aNAwI4W65TJ/OVeuXlJCM1HdWaRSFCUEm+/573KCkhSS1bxpBmwQIST79+tFiffELSOn2aFm7TJv524EBaqfHj2QZ+fgaRmO70q4QyN7f2FWymh9GonERODn/XtCnd0b172XbqVGmlsFlZVO6oKGNVYkwMr1GHzTz1FAnm22/53K1bs48++4zPa2dnroy7dxvDjP378z6//MJ7R0aaL3BS2fjAQCMES001pvKeOkWlPHaMdUpNNTZy+fVXktWiRVTsxYuZrc/IYDjp48Nnvvde1jslxdhAROWJ1DmU6rwFtS0cQHn69VcSQUwMR1+eeILPbzpd+N57Wf/585kHiY/n7376iTmxVavYHklJ/Oz9942pz82b05h0706ZSE/ndO+ahpKtRaMmAx8fY494d3c2+qhRjGHPnKHXcP48hWbdOlqnsDBafzXu6+HB9+qcABVTC0HFMh33PXuWCpycTKWIjTVi0LIyw5OozgJVRuvWzAgPGmQcchISQsHr2ZOvnZ2pUF27Goezengwhp8yxZjxqKC2RfP2Np/NqLbgUv8twdRzUC7u4MF8ff48FadvX7Zvp07GTMayMipJcLCxjDgjg++HDuU1W7awzv/4BwnmjTf43Y4dVMyHHiLhzZ3Lpdvffkurf/w4w6XRo+n19OxJMtywgaRi+izqRCW1Tb7aW/DqVbr3Xl5sTy8vY+WqCn0iI42DRyIj+awqFzF1qtE2arm0Ov4uPZ0K2b07CcI0cSiEQb6bN7NMtTqzXz8qvb+/kbg1XZE6eTLDYCcn3m/TJuZvTp9miLFxIz/Lz2f9Zs2ih9i+Pb3CgQON1bYTJtBw1heNeg9EX18KUno6rWfHjkxe+fqStQ8epNvp6WkMA2ZmUoCUKwxQKFaupCCq5B1gLN3t1o3vjx5lY1+7RlbeuNE4zj0lxZgjUHnRSnVo08Y4JLNTJyr6889zbHjYMHb20aNUhMpKbwkpKbzuiy/oJaWnU3hdXChI4eEsr7qlwpY8mpAQWqfLl83DphYtzKd0qzX6AIXv3Dm2wT33kOwefJDP1aEDE6iRkRT8Pn048qBi/MRE9lnnzoZ30awZE7zz59Nql5QwT1Q5wbZlC+Nw01DNw4NtsmsX+79DB36WnEx5SU42CEUND3fuzLAlMpLDzCNH0jAkJfH7w4dJaJcuMRwoLKTC79pFmZGSZa1ZYyz02rSJ9Rg2jMRhKaGfl2futQ0YQPIMD2ce4ssvSS6pqZS16dNJDGPGsH927KB89ulDcrB2i3Vr0ag9A7VD0V//Smtx9Cg7YtIk5gRKS+lSDRxIIVHLVivP4Vb7+6nZbQq9ehnxphot6NiRVnfAAKPjk5PZyarc7OyqU0ktQSmnlIZbfvEirV+7drTG6vQha6A29gwMpMBfuWIkIZUnVHn36MqoTAh+frTI/v4sz9KCmNJS852i1WrF2FgSkpMTlfDRR6nsa9aQPO67j/3y7LN85u3bSdZBQXTn8/ONRVMREdy1WQhj2bQp4uPZd5UXAglBb3DTJmOoU+0KdfUqXyuMGWPkLkaPpiuflcXna9OGXooQtLjNm5NcVq5k25w/T7IrKWGYVFxMGVD7aZSVVSVQU5SXG3KTkUHPKCmJu3Tt2sXw4OxZ5qjUSVkdOxqndqlNVGbPJonamgiARkwGpaVsmG7dKLBqdVfHjnSP27alqx4bawjIpEnsYDXeDpDZ1Rz5yhbdVGmio+kWOjgYLntmprHPoZMTBebs2aqkUhNUcqp5cz7H+fOG4k6YYOyIVBM8PIx9/8+fpyt/4QIFMzPTXEFqymc4OrL+pucNODrSa2nRwjj+rTLUYSMKLi70nkaPZrvFxfFZfvqJfTR2LCfUfPAByfy3v6XL/Ic/8LlXruQz+fgAzs6pOHBgCLZvT8NLL1HQMzPNp5UnJNArsnRMfPPmvNeTTxqLxNTio8pt4+FBpVW7P6sp5ps300NR8zm6d+cz//ijsU6kSxfK3n33sb1Gjzbk6Ztv+Pp3v1PTg83Xx0RG0ttQ2+3FxVFWDxygtxoTYxz8O2UKy+7UyZjOPWmScfZkQID50nlbotGSgYMDLb6zMxtWzQkYPNiYWz5rFoVEjel7erIRv/mGruDx47TEI0bwe9MNKhSCghgPnzvHuCshgb/t1InC7exsbNsdEUFLb22YoOq7bBk7dvBgY9fmuDgKlulin+oQEMBk1P79tLDBwYYFVweoWINWrSjgXSutKXVzo7Xq0YPKXRnXrpnfo0sXKlD79hRUHlFPlzk7m2HCSy/xu337OBEpKYn3/c1vSBynTtH1PXbsLeTn70Vs7Hx8+SUTlzt3GkvFDx1i31R30EhAgHG82YgRbCeFzMyqlnrcOHoRAJ8pIIDezIkTrKvaA7JHDxJNcDBDqebNjSneSUlUXD8/9sm6ddyAF2A9tmwxX3fh4UECPHiQSc/UVHopS5fScLm6sm1efJHh1nPPcWhVjaRt2mQcE3jkCGXgdqBR5wwUQkLInkFBRpKtqIj5hOHDKTy9e5Plx42jlYiJocIoIgAsu/edO1Pps7PpzirLWFREgR0+3OhY5UlYM+1ZQQiGGAEBhrt37BjjPhcXkkpWVs2bU6gt165epQAWFdFTUTF+SgrJJi6u6vFfpmjThlOYLa2taNqUmfz33jMfr1fDtqZo1YpKMHMm30dGkiDWrOGfry/bNTiYBLt9OwkxK4v3bt8eyM93xfXrRpLiwIFFOHCAh8M+/ngBzp3jvadNM0aAqnum55/na39/ErtaGaqWEptCTSb69FPDa8jPZz2bNqVS9utHBV2zhq/37mVSUg1vqo1qevWi0k6aZBiHli353bRp9LR27GAd+valhzh9Og2UmmTVsydzKQ8/XPUEqm7dKGenTxthWXa2bTYysYRG6xlUhtq6C2DjKXZ/8UUm6dSmG97e7HDTWWEKlpY3d+9O6+PrS7dNrTs4dsxY3WeKtm3NZ/JZg1GjjI7u25e5CfUsisyqw6FDFNbAQGOmoBpeUhudxMUZE2kqW31TeHpWn08IDDSE9Z//NOq3d2/VZcRC8HlM29LNjd6AOjFo0yYmYNXBILNn869XLxLFo48mIyDAmEoLuMHN7SG0bXsWx46RAAMD+fv166vfUdnOzpyYR4wwDtwFaBDi4mi9161jvN+iBes+aRKz+q+/znAtP58KmpHBPITa2DYjg209fDivyc1lnmXxYpLwE0+Yt01uLkdDDh1i+YWFHEp+911jsdTcuQyfFi9mrsXSUXSTJtHLCAuz/ii5+uCu8AwAxovr17NxXVwYKrz3Hhv/6afJ9IMHk5Ufe4yCpLbIrglubuwsKSlwTZoYJyj36lVVESIijKW9t4I2bcjssbGMD9XSZEunHEdH87uBAyk04eG8PiyMdSwooOC1aMH2GD689vs//XT137VoQQsfE0NCePZZY8SiMu67r+pnPXow9vXwqHoupKMjs+UREVROT09/eHh44vLlQtjbu6CsrBDl5Z7IzPS7Od26ZUtaaH9/Yy6Ammxj2lalpfSaLl/mn5Scc9K+PZWsSxcOIZoSYdu2xkiCgwOJ6rXX6Al8/z1JYcwYkmFoKBPNDg4kk7Q0Gp38fHpTysOIjeV1ERH87cKF7Ldt2/hbDw/Oh8jK4rNYyoGYYs8eyqbpBrrWTnq7Fdw1ZODry0ZWm2S88gqZ/oEH2IHr1pHhn3zSsO4rVtBtdnauvhE3baKwlZUZB3OqEQZ1cq4phDDc41tF69aMg9V+eWPHcpaaqYU4cIAC1KsXhSw315jcYnpG3/HjXNa7fr11Y81quNUShgxhGDF4MGPZP/+5+r0LLRGPEMzj7NhhbLB64AC9rKIiY1Zfz55Uxv370/Hss3Mwe/ZsLFq0GPv3pyIvj4rdpAmVRwgqhNrIw8WFZNy0KQlfbcnu48N+DAsjYZ4/z+9ND5k1p1uRgQAADnhJREFUhZ8fw0C1AW3LlkwAvvkm675+PT2d1avp7r/8MudJHD3KEG/PHj6riwu9gOJikrWaNTh2LMsH6DlMmUJiOXeOdap8hF9lZGWRRF54wXx+iLX5oVvBXUEGqjHatqW7lJNDa+XkxFg0Kood8fnn5nPrp0why8+cyU6qvIY+PZ0JQzUmnpZGsjl0iDkIU0Y2hTUWuCaMHs1k0tq1TBi5uDDG/PJLWrFLl0gYkZFUhB076P5WJqb27UkGpuHNrcDLi6MMdnbGTs179zLZ1aoVY2d1DFynThRIlXeQkqR27RqTo05ODDFeeokKFxREsj51igm6J57gMOSaNcCePWtu5iP+/e+FAIy9ExYsYPjn5UVlbNaMz9u6Nfs+N5f3PXGCROPjQ+JScxPCwxkumE70qQw1/TwggJ5GWBjj+UWLjGSktzc9zIIChgr330+ZGjSIVn7sWLaLSiJmZJAMnJ0ZBvTvT5LdtIm/qW2/DYVFi0gapgnQ3r1vX74AAIS05TGuViIqKkoequ4ki0rYsIH/LR1HDpCRT59mRyYmcvjqm2+M6aPx8YxjvbzYiX/4g/Hbf/6T3kZqKi3RhAlUhKgolllfpa8JRUWcZvrAA0bSLyeHf6bx4WefkchmzDD//YkTTOK1bEmBV2vmbwVnz9Ktfv1147PduzkFWR1kAxh7DFy9SpJSnzdpwvZViioEFdXentd98AGV09nZOCX7++9JPNUhJ4fnNhw/TsW3t6enVFzMZ+7QgWTj789yCwpISCp5qJYVd+xoPm27MgoLqdSPPGLkpV59lb8dO5YhzdmzVMQ1a5jL6NmT+Zvz53mfnBw+04UL7KtJk6j099xDYjt/nrJlzaxVgDKYmgo884x119cFQojDUkqLKdlG7xl4eNR8/pyTk7F2PySE+YKnn+awjVr5uHkzGd50L4GNG9nx7dpRMdUW5GVlVITKG6fYGs7OjGNXr6YgBgTQwpiOduzcSe/lySer/r5bN4MoVXhzq1Cz7UwxdCit5Jo1Rn7Aw6P2g1AU3NzodX3yCYW6a1cqxfHjVOjevWv+vacnc0LLlzOH4ebGsXp/fyr9yZP0iLy8OGrh6Egy6tOH/bpnDz2LGzc4sWfoUGMbO4AejbMzPRwfH4YHAwaQANUBqJ98QnI5d47lqZ2boqPZNl26kIQKCihrjzxCAzJoEI3QypUM80wXWdWGY8dY/ttvW/8bW6HRjyb4+lo+/aY6qEkbDzxghAzjxtGKKpd03z4K1hNPsPNHjqT7deoUla9Xr9ubqFEID2dd//1v8zUSUjKHEBPDMWdLdVETqY4fr/l4OGvg7Fx1mzSAw10xMRTQumL1auYf/vpXY4SjTRtaT9NDa2qCEAyjHn+cBPCf/xgHsTo48LWzM8f9T54kwa9cSWUqKuI8Ay8vWu39+zkfJTqabaaSr6mpLNvHhwYkO5teYXIyvZv8fMpfVha9DB8fDuX6+VFWHBwYOnz6KTeSdXdnHcrK6M2pmafW4MoVJhzffLP60OZ2otF7Bi1b1j1p8uijjLXGjWPiZ9Qo5g127uSU15wcTn5ZvZrXAkyaPfEEhS842PbPUR0iImidNm40TnjKyaGXo2bsVYdWragEte28ZA0sDUkKwfUCr7zCsXxr2qWwkIuUIiI4lFYZeXm1b1pbGeHhPML+T38i8b34IkcCtm8nGXh6Gvs7pKYyhGrWjG75hQvG8veuXams16+zvYuLjd2yy8tJUmppvJOTWqyUiitXZqJTp29x7Zofrl4lwXToQC9kyhT+ftcuGpOICJJDXXHuHL2Bf/2r5t29bycafc5ASuu2GrOEggIKcno6k4/cN45unNpDQHkLKSkUkDlz6n4fWyE7m0OMleffV4f8fMa01m64UhPUCT2WUFbGODoigqsPq7vm448p1K+/Xv0c/R076I7XtKlsTVi2zNigpkULKvS335IA1L6VarPSa9dICrm5/O7yZSqyGqVQ08JLS42krZMTQwA7OzWL9Vns3/8fBAU9g7CwTzBgAI1GWRlnA5aWUjaHDrV+Ilpl7N7NZPLf/lb9TlW2Qk05g0ZPBrZAYSFdxOxsCnx6OoVi2DAKyYkTFILqprxqEJs2UZn9/Y0psRkZ9E6KiozNW283srIYfrRoQcusXHE1M/DHH1kfFQr4+pJgs7OZDLxxgx6nuztDCTs7ehxOTpSNM2eARYtcUVZW9fAJR0cXrFxZgIAA4+i9W0VKCsMCX19jFuXtxm0jAyHEdADzAHQB0FtKaZWGNzQZWEJxMYfPvL0pTLdr8cd/Iy5fZq7F3p7t17//rVv6+uDkSSZRy8tpxZs2ZVLP398Ir7KymPc4eZKegtqf4fJlxuhOTrwmK4ufu7jwmZydU3H8+Mu4dGktSkvz4eDghk6dpuCeexbA3d3YY9DJiR6Fn5+xPZ6PT1VXX0qGJydPkoASEni/hx+2ftWqLXA7RxNOAJgK4D/1LKfB4eR0e4cO/5sREFD/iVe2QLdu/EtLo0dw5Yqx8Mvbm5ZfnampToG+epWeQGgoP2ve3NjjwHxZsD/mzvXE4sWFcHFxQXFxIQYP9sRHH5lvNlpQQAU/f57DiDk59EDKyph/cHLiffPyjC3P3N0t79dwp1EvMpBSxgOAaIjUu4ZGNfDz4zR1gBZYbYV29apx5mSHDiSIwEDrE3Tp6emYM4czJBcvXoxUC2Pcrq6cl2JpMVVeHsOVsjISz51KDFoLm+QMhBC7AbxcU5hQ6Xi1yPN1Xe2joaFRb9QrTBBCbAdgaSP2N6SU66ythJRyMYDFAHMG1v5OQ0OjYVArGUgpR9R2jYaGxt2PRj8DUUNDo2FQLzIQQkwRQqQA6AdgkxDiB9tUS0NDo6FR39GE7wF8b6O6aGho3EHoMEFDQwOAJgMNDY0KaDLQ0NAAoMlAQ0OjApoMNDQ0AGgy0NDQqIAmAw0NDQCaDDQ0NCqgyUBDQwOAJgMNDY0KaDLQ0NAAoMlAQ0OjApoMNDQ0AGgy0NDQqIAmAw0NDQCaDDQ0NCqgyUBDQwOAJgMNDY0K1HcPxA+EEKeFEMeFEN8LIZrZqmIaGhoNi/p6BtsAdJdShgFIBPBa/aukoaFxJ1AvMpBS/iilLK14Gw2gVf2rpKGhcSdgy5zBkwC22LA8DQ2NBoRNjlcTQrwBoBTA1zWUY3rW4i1VVkND4/ah3serCSEeAzABwD2yhlNc9VmLGhqNG/U6REUIMQbAHwEMkVLm26ZKGhoadwL1zRl8DMADwDYhxFEhxL9tUCcNDY07gPoer9bRVhXR0NC4s9AzEDU0NABoMtDQ0KiAJgMNDQ0Amgw0NDQqoMlAQ0MDgCYDDQ2NCmgy0NDQAKDJQENDowKaDDQ0NABoMtDQ0KiAJgMNDQ0Amgw0NDQqoMlAQ0MDgCYDDQ2NCmgy0NDQAKDJQENDowKaDDQ0NABoMtDQ0KhAfY9Xe6viaLWjQogfhRABtqqYhoZGw6K+nsEHUsowKWUEgI0A/myDOmloaNwB1Pd4tRyTt00A6PMQNDTuUtRrd2QAEEK8A+BRANkAhtW7RhoaGncEooZDkHiBFcerVVz3GgAXKeX/VFPOzePVAHQCkGBF/XwAZFpx3Z1EY69jY68f0Pjr2NjrB1hfxzZSyhaWvqiVDKyFEKINgE1Syu42KZBlHpJSRtmqvNuBxl7Hxl4/oPHXsbHXD7BNHes7mhBs8nYigNP1KU9DQ+POob45g78KIToBKAdwHsCc+ldJQ0PjTqC+x6vdZ6uKVIPFt7l8W6Cx17Gx1w9o/HVs7PUDbFBHm+UMNDQ07m7o6cgaGhoAGgkZCCHGCCEShBBnhBCvWvjeWQjxbcX3vwgh2jay+r0khDhVMTV7R8XISoOitjqaXDdNCCGFEA2eHbemjkKI+yva8qQQ4pvGVD8hRGshxC4hRGxFX49r4Pp9LoS4IoQ4Uc33Qgjxz4r6HxdC9KzTDaSUd/QPgD2AJADtATgBOAaga6VrngXw74rXMwF828jqNwyAW8XruQ1ZP2vrWHGdB4CfAUQDiGpsdQQQDCAWgFfFe99GVr/FAOZWvO4K4FwDt+FgAD0BnKjm+3EAtgAQAPoC+KUu5TcGz6A3gDNSymQpZTGAFQAmVbpmEoClFa9XAbhHCCEaS/2klLuklPkVb6MBtGqgulldxwq8BeB9AIUNWbkKWFPHpwEslFJeBwAp5ZVGVj8JwLPidVMAlxuwfpBS/gzgWg2XTAKwTBLRAJoJIfytLb8xkEEggIsm71MqPrN4jZSyFJz67N0gtbOufqaYBbJzQ6LWOgohegAIklJubMiKmcCadgwBECKE2CeEiBZCjGmw2llXv3kAHhZCpADYDOCFhqma1airrJqh3msTbABLFr7yEIc119wuWH1vIcTDAKIADLmtNbJwawuf3ayjEMIOwN8BPN5QFbIAa9rRAQwVhoLe1R4hRHcpZdZtrhtgXf0eALBESvk3IUQ/AF9W1K/89lfPKtRLTxqDZ5ACIMjkfStUdb9uXiOEcABdtJrcJVvCmvpBCDECwBsAJkopixqobgq11dEDQHcAu4UQ58B4cn0DJxGt7ed1UsoSKeVZcP1KMBoG1tRvFoDvAEBKeQCAC7gmoLHAKlmtFg2ZAKkm6eEAIBlAOxiJm26VrnkO5gnE7xpZ/XqAyafgxtqGla7fjYZPIFrTjmMALK147QO6vN6NqH5bADxe8bpLhaKJBm7Htqg+gTge5gnEg3UquyEfpIYHHAcgsUKh3qj4bD5oZQEy8EoAZwAcBNC+kdVvO4B0AEcr/tY3tjasdG2Dk4GV7SgAfAjgFIA4ADMbWf26AthXQRRHAYxq4PotB5AKoAT0AmaBSwDmmLTfwor6x9W1j/UMRA0NDQCNI2egoaHRCKDJQENDA4AmAw0NjQpoMtDQ0ACgyUBDQ6MCmgw0NDQAaDLQ0NCogCYDDQ0NAMD/B5qez9TryUwLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "with torch.no_grad():\n",
    "    # Initialize plot\n",
    "    f, ax = plt.subplots(1, 1, figsize=(4, 3))\n",
    "    \n",
    "    # Plot training data as black stars\n",
    "    ax.plot(train_x.numpy(), train_y.numpy(), 'k*', zorder=10)\n",
    "    \n",
    "    for i in range(min(num_samples, 25)):\n",
    "        # Plot predictive means as blue line\n",
    "        ax.plot(test_x.numpy(), output.mean[i].detach().numpy(), 'b', linewidth=0.3)\n",
    "        \n",
    "    # Shade between the lower and upper confidence bounds\n",
    "    # ax.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)\n",
    "    ax.set_ylim([-3, 3])\n",
    "    ax.legend(['Observed Data', 'Sampled Means'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulate Loading Model from Disk\n",
    "\n",
    "Loading a fully Bayesian model from disk is slightly different from loading a standard model because the process of sampling changes the shapes of the model's parameters. To account for this, you'll need to call `load_strict_shapes(False)` on the model before loading the state dict. In the cell below, we demonstrate this by recreating the model and loading from the state dict.\n",
    "\n",
    "Note that without the `load_strict_shapes` call, this would fail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "state_dict = model.state_dict()\n",
    "model = ExactGPModel(train_x, train_y, likelihood)\n",
    "\n",
    "# Load parameters without standard shape checking.\n",
    "model.load_strict_shapes(False)\n",
    "\n",
    "model.load_state_dict(state_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
